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Abstract 

A  review  of  progress  in  large  eddy  simulation  of  complex  re3w:ting  flows  is  presented.  Topics 
include:  the  governing  equations  for  low  Mach  number  combustion,  assumed  PDF  subgrid-scale 
models,  subgrid-scale  modeling  for  the  variance  and  dissipation  rate  of  a  conserved  scalar,  inflow 
and  exit  boundary  conditions  for  confined  swirling  flows,  simulation  results  for  isothermal  swirling 
flow  with  experimental  validation,  and  results  from  a  preliminary  reacting  flow  simulation. 

1.  Introduction 

The  coaxial  jet  combustor  is  an  idealized  gas  turbine  combustor  geometry  for  which  a 
number  of  experiments  have  been  performed  imder  a  variety  of  conditions,  both  reacting 
and  non-reacting.  Figure  1  shows  a  schematic  of  a  coaxial  jet  combustor  in  which  gaseous 
fuel  is  supplied  to  the  central  pipe,  and  swirhng  air  is  supphed  to  the  surrounding  annulus. 
The  two  streams  begin  mixing  after  a  sudden  expansion  and  then  react  in  a  swirl-stabilized 
diflFusion  flame. 

Because  of  the  availability  of  experiments,  this  configuration  was  chosen  as  the  context 
for  the  development  of  large  eddy  simulation  (les)  based  prediction  methods  for  complex 
reacting  flows.  It  is  expected  that  les  will  become  an  important  tool  for  analysis  of 
flows  that  are  not  easily  modeled  using  Reynolds  average  approaches  or  in  which  unsteady 
dynamics  of  the  large  scales  are  important  to  understanding  or  predicting  flow  behavior. 
Turbulent  reacting  flows  appear  to  satisfy  both  of  these  conditions  and  thus  represent  a 
promising  and  technologically  important  application  for  LES. 


Fig.  1  Schematic  of  a  coaxial  jet  combustor. 

Johnson  and  Bennett^  conducted  experiments  of  non-reacting  confined  coaxial  jets  us¬ 
ing  water  with  Reynolds  number  80,300  based  on  overall  mass  flow  rate  and  jet  diameter. 
Roback  and  Johnson,^  using  the  identical  configuration,  later  extended  the  measurements 
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to  include  swirl  in  the  annular  stream  with  a  swirl  number  of  approximately  0.41.  Sommer- 
feld  et  al.3  performed  experiments  in  a  slightly  different  configuration  with  particle-laden 
air  (Reynolds  number  52,400  and  swirl  number  0.47,  based  on  overall  jet  properties).  Owen 
et  ah'*  performed  a  series  of  reacting  experiments  with  methane  and  air  under  a  range  of 
flow  conditions. 

AkselvoU  and  Moin®’®  completed  a  large  eddy  simulation  study  of  the  non-swirling, 
non-reacting  coaxial  jet.  They  included  a  passive  scalar  in  their  calculations  to  investigate 
mixing  between  the  central  and  annular  jets  and  vahdated  their  results  against  the  exper- 
inient  of  Johnson  and  Bennett.^  The  present  study  focuses  on  the  effects  of  swirl  and  heat 
release  and  on  subgrid  models  for  turbulent  combustion. 

The  topics  and  issues  discussed  in  this  paper  include:  the  governing  equations  for  low 
Mach  number  combustion  (Sec.  2),  subgrid-scale  modeling  for  these  equations  including 
proposed  dynanuc  models  for  subgr id-scale  variance  and  dissipation  rate  of  a  conserved 
scalar  (Sec.  3),  the  numerical  techniques  used  in  this  study  including  a  novel  method  for 
generating  swirhng  inflow  boimdary  conditions  (Sec.  4),  results  for  isothermal  swirUng  flow 
including  the  effect  of  exit  boimdary  conditions  on  confined  swirling  flows  (Sec.  5),  and 
the  effects  of  heat  release  as  found  in  a  preliminary  reacting  simulation  (Sec.  6). 


2.  Governing  Equations 

We  consider  the  variable  density  momentum  and  conserved  scalar  transport  equations: 


^  +  (pun)  =  -  Vp  -h  V  .  [2p(S  -  |I  V  .  u)] 
at 


(1) 


dp^ 

dt 


+  V  •  (pu0)  =  V  •  {paV4>) 


(2) 


where,  S  is  the  strain-rate  tensor,  I  is  the  identity  tensor,  0  is  any  transported  scalar, 
and  a  is  the  molecular  diffusivity  of  the  scalar.  The  low  Mach  number  approximation  is 
invoked  so  that  only  the  uniform  background  pressure  po  enters  into  the  ideal  gas  equation 
of  state:  _ 

^  (3) 


This  effectively  decouples  the  density  from  the  pressure  and  precludes  the  formation  and 
propagation  of  acoustic  waves.  The  low  Mach  number  approximation  also  allows  the  energy 
equation  to  be  written  in  conserved  scalar  form  (2).  Since  the  density  is  determined  by 
the  equation  of  state,  the  continuity  equation  becomes  a  constraint  on  the  velocity  field, 
with  the  time-derivative  of  density  as  a  source  term: 


V(pu)  =  -|  (4) 

This  constraint  is  enforced  by  the  pressure,  in  a  manner  analogous  to  the  enforcement  of 
the  incompressibility  constraint  for  isothermal  flows. 

With  the  additional  assumptions  of  fast  chemistry  and  unity  Lewis  number,  the  den¬ 
sity,  as  well  as  temperature  and  all  chemical  species,  can  be  directly  related  to  the  conserved 
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scalar  mixture  fraction  Z,  which  satisfies  (2).  The  equation  of  state  then  reduces  to  a  den¬ 
sity  function,  p  =  p{Z).  Figure  2  illustrates  a  specific  form  of  the  dependence.  The  flame 
is  located  at  the  stoichiometric  value  of  the  mixture  fraction,  Zgt,  where  temperature  is  a 
Tnaximnm  and  density  is  a  minimum  (cf.  Ref.  7). 

When  the  fast  chemistry  and  unity  Lewis  munber  assumptions  are  relaxed  and  more 
general  combustion  models  such  as  laminar  flamelets®  are  used,  the  density  function  can 
be  extended  to  include  more  independent  variables  such  as  the  scalar  dissipation  rate  x- 
This  more  general  equation  of  state  is  written  as: 

P  =  p{Z,  X,  •  •  •)  (5) 


Fig.  2  Typical  dependence  of  density  and  temperature  on  mixture  fraction  from  the  fast  chemistry 
assumption. 

The  determination  of  the  continuity  constraint  source  term  —  ^  presents  some  diffi¬ 
culty.  It  can  be  determined  in  principle  by  differentiating  the  equation  of  state  with  respect 
to  time.  But  this  is  impractical  because  the  functional  dependence  (5)  can  be  arbitrarily 
complex  and  may  include  contributions  from  subgrid-scale  models  (discussed  in  Sec.  3) 
and  chemical  reaction  models.  The  typical  approach  is  to  approximate  the  time-derivative 
of  density  by  finite  difference  in  time. 

3.  Subgrid-Scale  Modeling 

The  LES  formalism®  introduces  a  filtering  operation  that  is  applied  to  the  governing 
equations  and  flow  variables  to  remove  the  imresolved  small  scales.  As  a  result  of  filter¬ 
ing  the  nonlinear  terms,  vmclosed  subgrid-scale  terms  analogous  to  the  unclosed  terms  in 
Reynolds  averaged  approaches  appear  in  the  equations.  Grid-filtering  is  indicated  using  an 
overbar,  and  a  tilde  is  used  to  indicate  density  weighted  (Favre)  grid-filtering,  (f>  =  ptj>/'p, 
so  that  is  the  resolved  scalar  field  computed  directly  by  LES. 

Subgrid  momentum  and  scalar  transport  terms  that  appear  in  equations  (1)  and  (2) 
are  modeled  using  the  dynamic  approach  of  Moin  et  al.^®  The  other  major  modeling 
requirement  is  for  the  nonlinear  density  function  (5).  One  needs  to  obtain  the  filtered 
density, 

P  =  p{Z,X,--)  (6) 
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but  the  scalar  transport  equation  provides  only  the  filtered  mixture  fraction  Z ,  and 


because  of  the  nonlinearity. 

A  general  method  of  evaluating  a  nonlinear  expression  such  as  (6)  is  to  express  it 
in  terms  of  the  probability  density  function  (PDF)  of  the  independent  variables.  In  the 
present  discussion  we  restrict  our  attention  to  the  case  of  fast  chemistry,  which  implies 
that  p  is  only  a  function  of  Z.  If  the  PDF  of  the  subgrid  mixture  fraction  fluctuation  P{Z) 
were  known,  the  filtered  density  could  be  evaluated  using 

p^-^  =  J  p(Z)P{Z)dZ  (7) 

The  PDF  can  be  modeled  directly  by  assuming  a  particular  analytical  form  for  it.  Prankel  et 
al.^^  and  Cook  and  Riley^^  proposed  using  the  Beta  family  of  distributions,  parameterized 
by  the  mean  and  variance.  Two  examples  of  Beta  distributions  are  plotted  in  Fig.  3. 
When  the  subgrid-scale  is  mostly  mixed  and  the  variance  is  small,  the  pdf  is  similar  to  a 
Gaussian  distribution.  But  when  the  variance  is  large,  the  distribution  becomes  bimodal, 
representing  a  subgrid-scale  state  that  is  mostly  unmixed.  Although  its  usefulness  in 
Reynolds  average  modeUng  has  been  questioned,  the  assumed  Beta  pdf  was  recently  shown 
to  accurately  represent  the  subgrid-scale  mixture  fraction  in  LES  by  Jimenez  et  al.^®,  who 
performed  a  priori  tests  in  a  turbulent  mixing  layer. 


Fig.  3  Typical  Beta  distributions. 

Once  particular  values  of  the  mean  and  veiriance  are  specified,  the  corresponding  Beta 
pdf  is  known  and  can  be  substituted  into  (7)  emd  integrated  to  obtain  the  filtered  density. 
In  this  way,  the  density  can  be  precomputed  and  stored  into  a  lookup  table,  to  be  retrieved 
during  the  simulation  as  a  function  of  the  resolved  mixture  fraction  field  Z  (the  mean), 
and  the  subgrid-scale  mixture  fraction  vauriance  Z""^, 


p  =  /(Z,  Z"2) 
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Fig.  4  A  lookup  table  for  density  as  a  function  of  Z  and  for  fast  chemistry  and  the  Beta 
PDF. 

An  example  of  a  lookup  table  is  shown  in  Fig.  4.  The  curve  for  zero  variance  is  the  same 
as  the  original  density  function  p[Z)  in  Fig.  2.  As  the  variance  is  increased,  the  subgrid 
field  becomes  increasingly  unmixed  and  therefore  unburned,  leading  to  higher  densities. 
The  mean  Z  is  known  firom  the  scalar  transport  equation  (2),  but  the  variance  requires 
additional  modeling.  A  model  for  the  subgrid  variance  of  a  conserved  scalar  is  discussed 
next. 

Subgrid-Scale  Variance 

To  predict  the  subgrid-scale  variance  of  a  conserved  scalar.  Cook  and  Riley^^  proposed 
a  scale  similarity  model,  which  they  showed  through  a-priori  tests  in  homogeneous  isotropic 
turbulence,  yielded  accurate  predictions.  The  model  has  since  been  used  in  an  LES  of 
nonpremixed  combustion  in  homogeneous  isotropic  turbulence^^  and  a  priori  tested  in  a 
turbulent  miYing  layer. The  major  drawback  of  this  approach  is  that  it  requires  input 
from  the  user  in  the  form  of  a  model  coefficient.  Furthermore,  there  is  no  reason  to  expect 
that  a  ‘universal’  value  for  the  model  coefficient  exists,  except  within  a  well  developed 
inertial  subrange.  In  general,  the  coefficient  could  vary  with  flow  type,  characteristics 
of  the  grid  and  test  filters,  Reynolds  and  Peclet  numbers,  etc.  Recently,  some  attempt 
has  been  made  to  determine  the  dependence  of  the  coefficient  on  such  parameters.^®  This 
Umitation  can  be  overcome,  in  principle,  by  using  the  Dynamic  Procedure^®’^®"'®  combined 
with  an  appropriate  scaling  law,  in  which  model  coefficients  are  determined  automatically 
in  response  to  local  flow  conditions.  _ 

For  the  purpose  of  modeling  subgrid  Favre  variance  of  a  conserved  scalar  ([>  the 
following  scaling  law  has  been  proposed^®: 

p^  =  C'A2plvJl2  (8) 

The  local  filter  width.  A,  serves  as  the  length  scale  of  the  subgrid  turbulence.  Eq.  (8) 
is  analogous  to  Yoshizawa’s  expression^®  for  the  subgrid  kinetic  energy  and  can  be  found 
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Fig.  5  Compcirison  of  predicted  subgrid-scale  variance  in  a  meridional  plane  using  the  (a)  scale- 
similarity  model,  and  (b)  dynamic  model. 


through  consideration  of  the  transport  equation  for  the  subgrid-scale  variance,  combined 
with  the  local  equilibrium  assumption. 

We  now  apply  the  Dynamic  Procedure  to  obtain  the  model  coefficient  C  in  (8).  To 
simplify  notation,  we  introduce  the  density-weighted  test-filtered  scalar  field,  defined 

by:  _ 

(9) 

P 

A  hat  indicates  test  filtering  at  a  scale  larger  than  the  grid  filter,  A,  usually  taken  to  be 
twice  the  grid  filter  width.  We  also  define  the  Leonard  term, 

£  =  (10) 

and  the  model  term. 

Ad  =  ^  1 V -  A2"^|vJp  (11) 

Note  that  the  scale  similarity  model^^  can  be  expressed  using  this  notation  as 

=  (12) 

where  Cs  is  the  unknown  scale-similarity  coefficient.  Next,  using  the  identity. 


P02  _  ^  _  ^p(j)2  _  P02  j  ^  £ 


6 


and  flJisnTning  that,  while  the  model  coefficient  may  vary  in  space  and  time,  the  same 
coefficient  apphes  to  both  filter  levels,  one  obtains  C  =  CM.  Although  C  could  be 
determined  locally  using  C  =  C/M,  least-squares  averaging  should  be  performed  over  the 
statistically  homogeneous  directions  to  obtain  a  more  accurate  (and  stable)  value.  For 
fully  inhomogeneous  flows,  the  more  general  Dynamic  LocaUzation  Procedure  of  Ghosal 
et  al.^®  should  be  used.  With  statistical  averaging  over  homogeneous  directions  indicated 
by  angle  brackets,  the  final  expression  for  the  model  coefficient  is 


{CM) 

(M^) 


(14) 


The  dynamic  and  scale-similarity  models  have  been  tested  in  the  large  eddy  simulation 
(presented  in  Sec.  6)  of  the  Roback  and  Johnson^  configuration  with  heat  release.  In  Fig.  5, 
the  performance  of  the  dynamic  model  is  compared  qualitatively  to  the  predictions  of  the 
scale  similarity  model  of  Cook  and  Riley. A  value  of  Cg  0.5  yielded  values  of  roughly 
the  same  magnitude  as  the  dynamic  model  and  was  used  to  obtain  Fig.  5a.  The  two 
model  predictions  are  similar,  but  the  dynamic  model  is  seen  to  provide  a  sharper  subgrid 
variance  field  withjnore  fine  scale  structure  due  to  its  direct  dependence  on  the  resolved 
scalar  gradient  |V(^|.  The  scale-similarity  model  generally  predicted  significantly  higher 
levels  of  variance,  though  the  overall  level  depends  on  the  value  chosen  for  Cg. 


Scalar  Dissipation  Rate 

The  dissipation  rate  of  the  variance  of  a  conserved  scalar  is  needed  as  a  parameter 
in  flamelet  models  of  nonpremixed  combustion.^^  Even  though  an  appropriate  scaling 
law  for  the  subgrid-scale  dissipation  rate  can  be  derived,  the  Dynamic  Procedure  cannot 
be  applied  directly  to  this  problem  because  virtually  all  of  the  dissipation  occurs  at  the 
smallest  unresolved  scales.  Instead,  the  approach  adopted  here  is  to  determine  the  subgrid 
variance  production  rate  and  assume  that  the  subgrid  dissipation  and  production  rates  are 
in  equilibrium.  A  similar  model  vfas  previously  proposed  by  Girimaji  and  Zhou,^^  though 
they  obtained  it  starting  from  different  assmnptions  and  using  spectral  transfer  arguments. 

The  assumption  that  the  subgrid  scEiles  are  in  local  equilibrium  with  the  large  scales 
leads  to  a  simple  model  for  the  subgrid  scalar  dissipation  rate  XSGS-  In  equilibrium,  XSGS 

is  equal  to  the  local  subgrid  variance  production  rate,  —p{uj<f/')<f>j  (repeated  indices 
imply  summation  and  a  subscripted  comma  indicates  partial  differentiation),  which  can 
be  obtained  by  assuming  a  gradient  diffusion  model  for  the  turbulent  scalar  flux  and  using 
the  dynamic  model  of  Moin  et  al.^®  for  the  turbulent  scalar-diffusivity  at  =  CaA^|.S|. 
One  then  obtains  xsgs  =  p(^t  |V<^p.  If  the  resolved  dissipation  rate  is  included,  the  total 
filtered  dissipation  rate  becomes 

X  =  p{a  +  at)  I (15) 

where  a  is  the  molecular  scalar-diffusivity.  Note  that  with  this  model,  x  can  be  obtained 
at  virtually  no  additional  computational  cost,  since  at  would  have  already  been  computed 
for  modeling  the  turbulent  scalar  flux  in  the  transport  equation  for  the  resolved  scalar 
field. 
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Fig.  6  Typical  grid  discretization  in  cylindrical  coordinates. 

4.  Numerics 

The  governing  transport  equations  (1)  and  (2)  are  solved  in  cylindrical  coordinates. 
Fig.  6  shows  a  typical  discretization,  which  includes  grid  refinement  near  solid  surfaces. 
The  numerical  scheme  consists  of  conservative,  second-order  finite-volumes  on  staggered 
grids.  A  second  order,  semi-implicit  time  eidvancement  approach,  in  which  the  derivatives 
along  a  single  coordinate  direction  are  treated  implicitly,  is  used  to  remove  the  CFL  re¬ 
striction  of  near- wall  grid  refinement  and  the  convergence  of  the  radial  coordinate  lines 
near  the  centerline.  Details  are  provided  by  Akselvoll  and  Moin.®  The  equation  of  state 
(3)  is  conveniently  and  efficiently  evaluated  from  the  independent  variables  (Z,  etc.)  by 
precomputing  the  density  function  and  storing  it  into  a  lookup  table.  The  continuity  con¬ 
straint  (4)  is  enforced  at  each  time  step  by  solving  a  Poisson  equation  for  the  pressure.  Fast 
Fourier  transforms  are  applied  in  the  azimuthal  direction  to  factor  the  three-dimensional 
Poisson  equation  into  a  set  of  two-dimensional  Helmholtz  equations.  The  Helmholtz  equa¬ 
tions  are  then  solved  with  either  a  line-relaxation  scheme  or  multigrid  scheme. 

Swirl  Generation 

Physically,  swirl  is  generated  by  tangential  jets  or  contoured  or  plane  vanes,  which  are 
difficult  to  simulate  numerically.  As  modern  simulation  and  experimental  methods  allow 
for  more  detailed  characterization  of  flow  behavior,  it  has  become  increasingly  important  to 
provide  well-defined,  reproducible  inlet  conditions.  This  problem  has  not  been  previously 
addressed  for  swirhng  flows. 

Computation  of  turbulent  swirling  flows  has  been  carried  out  mostly  in  the  Reynolds 
averaged  turbulent  flow  prediction  context.  There,  the  inflow  condition  is  usually  deter¬ 
mined  by  specifying  either  experimental  data  or  else  simple  algebraic  mean  profiles.  Further 
discussion  is  provided  in  a  review  article  by  Sloan,  Smith,  and  Smoot. The  present  study 
appears  to  be  the  first  attempt  at  large  eddy  simulation  of  tmbulent  swirling  flows.  The 
azimuthal  body  force  technique  described  here  represents  a  means  of  predicting,  rather 
than  prescribing,  swirling  inflow  boundary  conditions. 


8 


x/R  AXIAL  VELOCITY  SWIRL  VELOCITY 


Fig.  7  Comparison  of  predicted  radial  profiles  of  axial  and  azimuthal  velocity  with  experimental 

data  of  Sommerfeld  et  al.^: - LES  226  x  163  x  128  grid,  corresponding  to  the  bottom  of  Fig.  9; 

. LES  corresponding  to  the  top  of  Fig.  9;  •  experiment. 


Equilibrium  swirling  flows  can  be  generated  numerically  by  solving  for  the  flow  driven 
by  fictitious  axial  and  azimuthal  body  forces  in  a  spatially  periodic  pipe,  from  which  data 
is  fed  as  inflow  boimdary  condition  into  the  main  computational  domain.  The  axial  body 
force  represents  the  mean  pressure  gradient  which  drives  the  physical  flow  and  has  long 
been  established  as  a  means  of  driving  spatially  periodic  pipe  and  channel  flows.  On  the 
other  hand  the  azimuthal  body  force,  used  to  drive  the  swirl  component,  is  not  physically 
producible  and  should  be  thought  of  as  existing  only  to  overcome  drag  from  the  walls.  The 
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r/R 

Fig.  8  Comparison  of  predicted  radial  profiles 
mental  data  of  Roback  and  Johnson  :  - LES 


r/R 


of  mixtine  fraction  (passive  scalar)  with  experi- 
250  X  125  X  64  grid;  •  experiment. 


resulting  flow  is  a  stationary  approximation  to  slowly  decaying  swirl. 

The  radial  profile  of  the  axial  body  force  /®(r)  is  taken  to  be  uniform  as  in  the 
Poiseuille  flow.  For  the  swirling  case  there  does  not  exist  a  corresponding  physical  mech¬ 
anism  for  driving  the  flow,  so  there  is  some  freedom  in  choosing  the  forcing  profile  fe{r). 
Tests  revealed  that  the  flow  is  rather  insensitive  to  the  forcing  distribution  used.  A  uni¬ 
form  body  force  appeared  to  be  adequate,  and  for  simphcity,  it  was  chosen  for  use  in 
all  subsequent  simulations.  Details  and  further  discussion  will  appear  in  a  forthcoming 
article.^"^ 


5.  Isothermal  Flow  with  Swirl 

The  equations  governing  isothermal  flow  are  Eqs.  (1),  (2),  and  (4)  with  p  =  const. 
We  consider  the  incompressible  flow  experiments  of  Roback  and  Johnson^  and  Sommerfeld 
et  al.^  Figures  7  and  8  present  a  comparison  of  the  computed  mean  velocity  and  mixture 
fraction  profiles  with  experimental  data.  Although  a  few  discrepancies  can  be  seen,  overall 
agreement  is  good.  In  particular,  note  the  agreement  at  the  axial  station  closest  to  the 
inlet  which  shows  the  effectiveness  of  the  above  swirl  generation  method. 

While  the  Reynolds  average  approach  has  been  fairly  successful  at  predicting  the  ve¬ 
locity  field  for  swirling  flows,  it  has  been  much  less  successful  at  predicting  scalar  mixing.^^ 
The  LES  approach,  on  the  other  hand,  can  give  very  accurate  predictions  as  evidenced  in 
Fig.  8,  probably  due  to  its  accurate  representation  of  the  large  scale  mixing. 

Exit  Boundary  Conditions 

The  simulations  have  shown  that  exit  boundary  conditions  can  have  a  tremendous  in¬ 
fluence  on  confined  swirling  flows.  The  presence  of  a  contraction  or  expansion  at  the  end  of 
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the  test  section  can  completely  change  the  character  of  the  upstream  central  recirculation 
zone.  This  has  been  observed  in  experiments^®  and  has  been  known  in  vortex  breakdown 
computations,  but  was  not  previously  reported  in  turbulence  simulations.  Early  simula¬ 
tions  using  the  convective  outflow  boundary  condition^^  produced  the  flow  depicted  at  the 
top  of  Fig.  9,  contrary  to  the  experimental  results.®  In  later  simulations,  the  geometry  was 
modified  by  adding  a  second  expansion  at  the  outflow  boundary  to  represent  the  effect 
of  a  stagnation  chamber  at  the  end  of  the  experimental  test  section.®  This  modification 
produced  a  very  different  flowfield,  shown  in  the  bottom  of  Fig  9,  that  compares  well  with 
the  experiment.® 


Mean  Streamfunction 


Fig.  9  Effect  of  downstream  boundary  conditions  on  a  confined  swirling  flow.  The  character  of 
the  central  recirculation  zone  is  completely  changed  by  the  addition  of  a  second  expansion  to  the 
geometry. 

Perhaps  the  best  way  to  understand  this  phenomenon  is  to  consider  the  behavior  of  the 
pressure  in  a  swirling  flow.  Whenever  there  is  swirl  or  rotation  in  a  flow,  a  radial  pressure 
gradient  must  exist  to  balance  the  centrifugal  force.  Without  an  exit  expansion  a  strong 
vortex  core  naturally  develops,  thereby  decreasing  the  pressure  along  the  centerline  and 
creating  a  favorable  axial  pressme  gradient  which  may  further  enhance  the  development  of 
the  vortex  core.  The  effect  of  the  exit  expansion  is  to  force  the  pressure  to  become  radially 
uniform  at  the  exit  plane,  causing  an  adverse  axial  pressme  gradient  along  the  centerline 
that  may  influence  the  development  of  the  flow  all  the  way  back  to  the  inlet. 

Escudier  and  Keller®®  studied  the  effect  of  exit  contraction  on  confined  swirling  flows 
but  did  not  consider  whether  exit  expansion  would  have  any  effect.  They  also  speculated 
that  reacting  flows  would  not  be  susceptible  to  the  exit  boundary  condition  because  of 
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substantially  reduced  density  and  increased  axial  velocity  due  to  heat  release.  Apparently, 
Reynolds  average  calculations  of  confined  swirling  flows  have  not  shown  such  sensitivity 
to  exit  boundary  conditions. 

6.  Reacting  Flow  with  Swirl  and  Heat  Release 

As  a  first  step  towards  simulation  and  modeling  of  real  reacting  systems,  the  capability 
of  the  variable  density  formulation  to  capture  effects  of  heat  release  was  tested  in  an 
idealized  combustion  system,  by  assuming  infinitely  fast  reaction  of  the  fuel  and  oxidizer 
streams  in  the  Roback  and  Johnson^  configuration.  The  present  results  are  preliminary  in 
that  they  have  not  been  validated  against  experimental  data. 


Fig.  10  Comparison  of  mixture  fraction  fields  for  (a)  non-swirling  isothermal  flow,  (b)  swirling 
isothermal  flow,  and  (c)  swirling  flow  with  heat  release  in  identical  configurations.  Black:  air, 
white:  fuel,  and  gray:  mixed. 
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The  combustion  parameters  for  this  case  were  chosen  to  reflect  typical  conditions  in  a 
practical  combustor.  The  density  ratio,  Pmax/Pmin,  was  taken  to  be  3,  which  is  consistent 
with  preheated  air  from  the  compressor,  and  the  equivalence  ratio  was  set  to  0.9.  Fig.  10 
compares  snapshots  of  the  mixture  fraction  for  this  case  and  two  previous  isothermal  cases, 
all  from  the  same  Johnson-Bennett-Roback^’^  configuration.  In  the  non-swirling  isothermal 
case,  the  fuel  and  oxidizer  streams  mix  very  gradually.  When  swirl  is  added,  the  mixing 
effectiveness  of  the  flow  increases  dramatically.  When  heat  release  is  added,  the  mixmg 
region  is  pushed  outward  as  the  flow  accelerates,  and  the  mixing  efficiency  is  substantially 
reduced.  The  annular  air  flows  around  the  central  mixing  zone  and  impinges  on  the  outer 
wall,  forming  a  wall  jet  which  acts  to  confine  the  core  flow. 

Fig.  11  shows  a  quantitative  comparison  of  the  mean  velocity  profiles  for  the  isothermal 
and  reacting  cases.  Near  the  expansion  point,  before  much  heat  release  has  occiirred, 
the  profiles  are  very  similar,  but  further  downstream,  one  can  see  the  axial  flow  being 
accelerated  by  the  heat  release.  On  the  other  hand,  the  swirl  component  of  velocity  is 
largely  unaffected  by  the  heat  release,  as  can  be  expected  from  the  fact  that  the  density 
variations  aire  axisymmetric  on  average. 

7.  Conclusion 

Large  eddy  simulation  of  a  coaxial  jet  combustor  with  swirl  was  performed.  The 
velocity  and  concentration  profiles  are  in  good  agreement  with  experiments.  It  was  shown 
that  confined  swirling  flows  can  be  very  sensitvie  to  downstream  boundary  conditions.  To 
account  for  the  effects  of  heat  release  a  low  Mach  number  variable  density  formulation  for 
LES  was  developed  and  implemented.  An  assumed  PDF  approach  was  used  to  model  the 
subgrid-scale  mixture  fraction  in  the  limit  of  fast  chemistry,  and  a  new  dynamic  model  was 
used  to  compute  the  variance  of  the  subgrid  mixtmre  fraction  as  input  to  the  assumed  PDF. 
Simulation  results  showed  that  heat  release  reduces  mixing  and,  as  expected,  accelerates 
the  flow  in  the  combustor. 


8.  Future  Directions 

Subgrid  models  to  account  for  finite-rate  chemistry  need  to  be  developed  and  validated 
against  combustion  experiments  such  as  Owen  et  ah'*  The  benchmark  experiments  of  a 
piloted  methane/air  jet  (Barlow  et  al..  Combustion  Research  Facility,  Sandia  National 
Laboratory)  will  be  simulated  to  provide  detailed  validation  of  the  chemistry  model.  An 
extension  of  the  low  Mach  number  equations  to  incorporate  low  frequency  acoustics  will  be 
used  to  capture  acoustic  combustion  instabilities.  The  resulting  simulation  capabilities  will 
be  used  to  develop  active  control  strategies  to  suppress  instability  and  improve  combustor 
performance.  It  is  hoped  that  this  effort  will  continue  under  support  from  the  AFOSR 
propulsion  program. 
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